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ABSTRACT 

We study relativistic unmagnetized collisionlcss shocks using unprecedentedly large particle-in-cell 
simulations of two-dimensional pair plasma. High energy particles accelerated by the shock are found 
to drive magnetic field evolution on a timescale > 10 4 plasma times. Progressively stronger magnetic 
fields are generated on larger scales in a growing region around the shock. Shock-generated magnetic 
fields and accelerated particles carry > 1% and > 10% of the downstream energy flux, respectively. 
Our results suggest limits on the magnetization of relativistic astrophysical flows. 
Subject headings: shock waves — magnetic fields — acceleration of particles — gamma rays: bursts 



Due to the low plasma densities, shock waves ob- 
served in a wide range of astronomical systems are col- 
lisionless, i.e., mediated by collective plasma instabil- 
ities rather than by binary particle collisions. Such 
shocks pla y a central role in, fo r exam ple, supernova 
remnants (iBlandford and Eichlerl Il987f) . jets of radio 
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hj|2003|), 7-ray bursts (GRBs, 



12005(1 . pulsar wind nebulae (PWN. lKirke t al. 2007)), and 



Piran 



the formation of large-sc ale structure in the Universe 
(|Loeb and Waxmanl l2000t l . It is widely accepted that 
particles accelerated to high energy in such shocks gen- 
erate the nonthermal radiation observed in a wide range 
of astrophysical sources and constitute the observed pop- 
ulation of cosmic rays. 

Despite intense research, collisionless shocks are still 
not understood from first principles. In particular, there 
is no self-consistent theory describing the acceleration of 
particles and the generation of magnetic field fluctua- 
tions, which in turn scatter particles and mediate their 
acceleration. Much of the research has focused on "mag- 
netized" shocks, where the upstream magnetic energy 
flux constitutes a significant fraction of the total energy 
flux. Here we focus on "unmagnetized" shocks, where 
the upstream magnetic energy flux is small. An ex- 
treme example of such shocks are the relativistic GRB 
afterglow shocks, where the magnetic fraction of the en- 
ergy flux, eb, increases from ~ 10~ 1 in the upstrea m to 
es — 0.01 — 0.1 in the downstream (| Waxmanl l2006t and 
references therein). Under such conditions, it is likely 
that the initial upstream magnetic field does not play 
a role in the de termination of the shock structure (e.g., 
lGruzinovl[2li0l . 

A near equipartition field, eb — 0.1, may be produced 
by electromagnetic { e. g., W ei bel) instabilities {e.g., 
Blandf ord and Eichler 19 8?| iGruzinov and Waxmanl 
19991: iMedvedev and Loeblll999( l. The coherence length 
of the generated field is expected to be comparable in 
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this case to the plasma skin-depth, l sc i — c/lu p , where 
Up is the plasma frequency and c is the speed of light. 
The main challenge associated with the downstream 
magnetic field is that while magnetic power on l sc i scales 
rapidly decays after the shock, observations imply that 
near eq uipartition fields must persist over lO 10 ^ down - 
stream (|Gruzinov and Waxmanl [T999t IGruzinovl 2001). 
This suggests that the m agnetic field develops similarly 
large coherence lengths ([Gruzinov and Waxmanl 119991 : 
IGruzinovl 120011 ). Such evolution could be driven by 
large scale currents carried by particles accelerated 
in the shock, possibly leading to a self-similar plasma 
configuration (Katz ct al. 2007). The mechanism for 
the generation of large-scale magnetic fields in shocks 
remains unknown. 

Particle-in-cell (PIC) simulations, in which the plasma 
is represented by macroparticles and Maxwell's equa- 
tions are solved on a grid, have been extensively 
used in recent years to study shocks. Such studies 
have quantified the generation of upst ream current fil- 
aments by pinching instabilities (e.g., [S ilva ct al. 120031: 
Frederiksen et al.ll2004lJaroschek et al.l20 05: Spitkovskv 
20051 l2008aHChang et al.ll2008l ). and resolved the forma- 



tion of shocks in two- and three - dimensiona l (2D and 3D) 
pair plasma (jSpitkovskvl 120051 : lKatol l2007t IChang et all 
2008) and in 2D ion-electron plasma (Spitkovsky 2008a). 
These simulations revealed rapid decay of magneti c 
fields downstream (jGruzinovl l200l IChang et all 12008( 1. 
leaving the question of field survival over scales 3> 
l s d open and trigger ing alternative suggestions for 
field generation (e.g., [Goodman and MacFadvenl 120071 : 
iMilosavljevic et al.lf2007t ISironi and Goodmanll2007( l. ~ 
In this Letter, we report new PIC shock simulations 
performed on unprecedentedly long length and time 
scales, (L/l s d) 2 (Tcu p ) ~ 4xl0 10 . These simulations show 
the growth of magnetic power on progressively longer 
scales driven by the accelerated particles, and impose 
lower limits on the efficiencies of particle acceleration and 
magnetization. Our results suggest that e ven the most 
exten sive simulations previously reported (jChang et al.l 
2008) were too small to capture significant particle ac- 
celeration and th e resulting magnetic field evolution (see 
iKatz et all 120071 . for a discussion of PIC simulation re- 
sults and limitations). Here we discuss only the main 
properties of shock evolution; the particle acceleration 



2 



mechanism is discussed separately (jSpitkovskvl [2008b) . 
and a detailed analysis of the simulations is deferred to 
a later publication. 

Simulation set-up. For simplicity, we focus here on 
strong, relativistic shocks in unmagnetized pair plasma. 
In order to reach long length and time scales, we resort 
to 2D, and comment below on expected differences with 
respect to 3D shock s. We use the ele ctromagnetic PIC 
code TPJSTA N-MP flSpitkovskvl l2005fl. a parallel version 
of TRISTAN (|Bunemanlll993[ ) heavily modified to min- 
imize noise and numerical instabilities. A rectangular 
simulation box is set up in the x — y plane, with periodic 
boundary conditions in the y-direction and a conducting 
wall at x wa u = 0. Cold, neutral plasma is continuously 



injected from x^ 



t < lOOOw- 1 , 



tions are sufficiently large to reveal the onset of particle 
acceleration and the slow evolution of shock properties 
(evident on ~ lOOOw" 1 timescales) driven by these en- 
ergetic particles. A small fraction of particles, acceler- 
ated to Lorentz factors 70 < 7 < Jmax by repeated 
scatterings near the shock, gradually builds up a flat 
(7 2 dn I dj ~ const.) power-law energy tail downstream, 
already contai ning a fraction e ac c > 10% of the energy 
at t = 10 4 w- 1 ()Spitkovskvll2008oh . Here we defined e acc 
as the ratio between the energy density of particles with 
7 > 570 behind the shock, and the far upstream kinetic 
energy density, such that for a thermal distribution of 
the particles in the downstream e ac .c(ln » 1) — 0.3%. 



ct in the — x direction, where t is 



the simulation time. Reflection off the wall then results 
in a shock propagating along +x. All parameters are 
measured in the downstream frame, in which the wall is 
at rest. 

Typical simulation parameters are: injected bulk 
Lorentz factor 70 = 15, thermal spread A70 = 10~ 4 , 
and N ppc — 8 particles per species per cell, with spatial 
and temporal resolutions Sx — l sd /W and St = 0.045a;~ 1 . 
Here uj 2 = 47r(n e + + n e - )q 2 /jo r n, where m and q are the 
particle mass and charge, and n is the upstream number 
density. Our largest simulation has ~ 2 x 10 10 particles 
and (L x /l sd ) x (L y /l sd ) x (Toj p ) = 6300 x 1024 x 6300 ~ 
4 x 10 10 , although smaller simulation boxes have been 
evolved for as long as 12600w~ . The results displayed 
below mostly refer to a simulation with L y = A02l sd , 
evolved for T = 12600w, 
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Short time evolution. At early times, 
we recover shock formation as reported previousl y 
(lSpitkovskvll2"005l : iChang et alJl20M lSpitkovskvl[2008ah : 
a transition layer of a few 10^ thickness propagat- 
ing upstream, in which the plasma isotropizes, ther- 
malizes and compresses. The simulated shock tran- 
sition agrees to within a few percent with (magnetic 
free) hydrodynamic jump conditions: a shock velocity 
v sh = c(r d -l)[(7 -l)/(7o + l)] 1/2 and density compres- 
sion ratio n d /n u = (T d +-f ~ 1 )/(T d - 1). Here, T d ~ 3/2 is 
the downstrea m adiabatic index, and upstream pressure 
was neglected (|Spitkovsky||2008ar ) . 

Upstream, the interaction between the unshocked flow 
and a counterstream running ahead of the shock leads to 
the formation of current filaments (in both 2D and 3D) 
parallel to the flow, surrounded by near-equipartition 
filamentary magnetic (in the fluid frame) structures. 
Behind the shock, near equipartition magnetic clumps 
form and are advected with the downstream flow in 
2D. At early times (where 3D simulations are possi- 
ble, t < lO 3 ^" 1 ), good agreement is found between 
these clumps and the 2D projection of extended mag- 
netic loops formed nearly perpendicular to the flow in 3D 
shocks. When averaged along the transverse direction, 
es = (i? 2 /87r)/[(7o — \)nmc 2 ] (where B is the magnetic 
field amplitude) peaks at ~ 7% near the shock transition 
layer and decays be low 0.1% within 1000^ downstream 
(|Chang et al.ll2008h . 

Long time evolution. The above description does not 
include the effects of high energy particles accelerated by 
the shock, negligible at early times. Our present simula- 




FlG. 1. — Plasma evolution within 1000i s[ j of the shock. Nor- 
malized transverse magnetic field sign(B)eg (color scale stretched 

in proportion to e^ 4 to highlight weak features) is shown at (a) 

early (t\ = 2250OJJT 1 ), and (b) late (£2 = 11925u;p 1 ) times. Here 
Ax = x — x sh is the distance from the shock, with x sh (dashed) de- 
fined as median density between far upstream and far downstream. 
Also shown are the transverse averages (at t\, dashed blue, and <2, 
solid red) of (c) electromagnetic energy normalized to the upstream 
kinetic energy, €em = [{B 2 + i? 2 )/87r]/[(7o — ljmnc 2 ] (with E the 
electric field amplitude, included because in the simulation frame 
the induced E ~ B upstream), (d) density normalized to the far 
upstream, and (e) particle momentum jI3 x (with (3 the velocity in 
c units) in the x-direction averaged over all particles (higher ( , y/3x)) 
and over downstream-headed particles only. 

The energetic particles running ahead of the shock sig- 
nificantly alter the properties of the counterstream and 
the resulting current filamentation and magnetization 
upstream. Figure [1] shows the resulting spatial distri- 
bution of magnetic fields at early vs. late times, as well 
as the density and momentum profiles. It reveals an 
increasing magnetization level, with fields generated on 
gradually larger scales and extending farther away from 
the shock, both upstream and downstream. As a result, 
the shock compression transition layer (defined, say, be- 
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tween 10% and 90% of full shock compression) widens, n 
and eb become more oscillatory with distance behind 
the shock, the shock slightly accelerates (by < 1%), 
and the final compression ratio slightly decreases (by 
< 4%). Due to the substantial energy carried by the 
accelerated particles running ahead of the shock, the 
average momentum is strongly modified far upstream, 
Ax = x — x s h > 1000i s d, although the incoming flow 
slows down considerably only at Ax < 1001 s d (figure le). 

Figure quantifies the evolution of downstream mag- 
netization. It shows the power spectrum Pk defined by 
eb = J Pkd log k (curves) , the average coherence length 
(A) and the energy fraction eb (filled circles) of magnetic 
fields in an l x — S00l s d long region trailing behind the 
shock. In order to avoid shot noise contamination 1 , eb is 
measured only for coherence lengths A = 2ir/k > 101 s d- 
During 10 3 < tui p < 10 4 , eb in this region grows by a 
factor of ~ 4 and reaches ~ 1%. As illustrated in Figure 
HI the typical size of upstream filaments increases sub- 
stantially in time, with typical thickness (10 — 15)l s d at 
t = lO 3 ^! 1 growing by a factor of 3 — 4 by t = lO 4 ^" 1 . 
Figure [2] shows a modest effect downstream, where the 
average magnetic coherence length grows by ~ 10% dur- 
ing this epoch. 




Fig. 2. — Magnetic 2D power spectrum in a sample down- 
stream region denned by —1000 < Ax/l a d < —200. The spec- 
trum is shown to gradually grow and possibly flatten with time 
[solid curves for tu) p ~ 1900 (blue), 4600 (green), and 12600 
(red)]. The integrated energy fraction eg and average scale (A) = 
1-k I expfe" 1 J log(fc) Pfe rflog k] of the magnetic field (filled circles) 
grow correspondingly. Suppressing particle acceleration (cooling all 
particles above 7 coo ; = 80, dashed line and square for tuj p = 5750) 
stops the magnetic evolution. Also shown are estimated shot-noise 
power (dotted) and a A = rj, ~ (2P^) _1 / 2 i s[ j curve (dash dotted) 
with rjj the Larmor radius of 7 = 70 particles, roughly separat- 
ing magnetized (large scale) and non-magnetized (small scale) bulk 
plasma regimes. 

In order to test the role played by high energy parti- 
cles in the evolution of the shock, we have performed a 
suite of simulations with artificially suppressed particle 
acceleration. In these runs, we introduced cooling, where 

Noise with power inversely proportional to N ppc , filtered on 
small scales, and tested not to distort our results. 



particles with 7 > "f coo i lose a random fraction of their 
nonthermal energy, with various choices of 7 CO o/ ■ Cooling 
is found to significantly slow down or completely stop 
shock evolution, leading to a fixed magnetization level 
and a steady-state magnetic power spectrum. Higher val- 
ues of "fcooi are found to produce larger eb and (A). For 
example, a steady state configuration with eb = 0.2% 
(at -1000 < Ax/hd < -200), obtained for -y cool = 80 
(70 = 15), is shown as a dashed line in Figure El 

The decay rate of magnetic fields advected downstream 
slows down as the shock evolves, as illustrated in Figure 
03 This is partly attributed to the increased fraction of 
power deposited in large scale fields, which are expected 
to survive farther downstream. Indeed, the inset of Fig- 
ure [3] shows that for A > 40Z s d, magnetic evolution is well 
fitted by exponential decay, Pk oc e _wt , with uj oc k 2 . 
Such behavior is expected, for example, in MHD mag- 
netic diffusion with scale-independent resistivity (note 
that Figure [2] suggests that the bulk plasma is indeed 
magnetized on these scales). Figure |3] thus indicates, for 
example, that power on A = 30^ (100l s d) scales survives 
more than 10 3 l s d (lO 4 ^) downstream. 




10" 2 IJk 10" 1 
10 -> I ^ , , J 

500 1A 1000 ,„ 1500 2000 

lAx-Ax J/1 J 

sd 

Fig. 3. — Decay of magnetic energy in three different slices ad- 
vected downstream at progressively later times. Each slice is l x = 
1000i s( j long, located with its upstream edge at Axo = — 500i s[ j 
at tou} p = 3600 (solid red), 5175 (dashed green) and 6750 (dash- 
dotted blue); all curves cutoff when the simulation terminates at 
tu) p = 8550. The decay rate is seen to slow down as the shock 
evolves. Power on scales A > 40i s[ j decays nearly exponentially in 
I Ax I, shown (dotted red) for touj p = 3600 andA/Z st j = 40, 50 and 
67 (fast to slow decay). In Chang ct al. (2008), magnetic decay 
was studied in a fixed downstream frame region (e.g., solid line), 
not taking into account magnetic field evolution due to accelerated 
particles. 

Inset: Spatial decay rate of the magnetic power spectrum 
(2/v sh )uj(k) = P~ 1 (-dP k /dAx) in a slice (with l x = 1400i sd 

and to = 4500a; p ) advected downstream at (t — to)u) p = 
0, 225, 450, ... 1800 (solid, dark to light curves). At large scales, 
the decay rate is well fitted by u> oc fc 2 (dashed), whereas at small 
scales ut decreases with time/distance from the shock. 

Convergence. Convergence tests were performed with 
respect to all simulation parameters (around their val- 
ues given above), with no qualitative changes to the re- 
sults. Figures Q] and [H indicate that the simulation box 
used is sufficiently large to avoid significant boundary ef- 
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fects. However, as ~ 40% (5%) of the magnetic power 
is already deposited in A > 50l s d (A > 100Z S( j) scales at 
t ~ 10 oC 1 , increasingly larger simulation boxes, in both 
longitudinal and transverse dimensions, will be required 
in order to properly resolve the shock and the growing 
coherent structures. 

Conclusions. Our analysis shows that collisionless 
shock configurations simulated previously may represent 
steady state configurations only as long as particle accel- 
eration remains insignificant. We find that a population 
of energetic particles is accelerated and drives the gener- 
ation of progressively stronger fields on gradually larger 
scales. Our simulations do not reach a steady state; 
rather, an increasing fraction of shock energy is trans- 
ferred to energetic particles and magnetic fields through- 
out the simulation time domain. 

Once stochastic acceleration and magnetization ensue, 
they are unlikely to diminish to a lower energy steady 
state. Hence, our results suggest lower limits to the 
efficiencies of magnetization and particle acceleration, 
£ b ^ 1% & t distances \Ax\ < D = 1000Z S( 2 downstream 
of the shock, and e acc > 10% with no significant cool- 
ing identified downstream. We find no evidence for the 
saturation of e acc , eg, D or j m ax, although the high en- 
ergy particles downstream are already sub-equipartition 
at t ~ lO 4 ^- 1 . 

Although our results are obtained for 2D pair plasma, 
we expect qualitatively similar shock evolution in 3D 
shocks and in electron-ion plasma. While the nature 
of upstream current filaments and downstream magnetic 
loops/clumps may depend on dimensionality, 2D and 
3D simulations are in good agreement at early times 
(Spitkovsky & Arons, in preparation). Also, shocks in 
ion-electron plasma were found to be similar to pair 
plasma shocks at earl y-time 2D simulati ons due to effi- 
cient electron heating (|Spitkovskvil200 8a). Some level of 



stochastic particle acceleration is inevitable in all cases, 
but the generalization of our results to 3D or to ion- 
electron plasma is yet to be tested at late times and in 
the presence of a high energy particle tail. 

The major role played by high energy particles in shock 
evolution, their flat spectrum and the apparently flatten- 
ing magnetic power spectrum are t rends consistent with 
a self-similar plasma configuration (jKatz et al.ll2007f ). al- 
though the simulated downstream scale growth is more 
modest than the A oc D self-similar scaling. At this stage, 
the simulations are not yet sufficiently advanced to vali- 
date or rule out self-similarity. 

In summary, we have shown that collisionless shocks 
in 2D pair plasma evolve on long, > lO 3 ^" 1 timescales, 
such that the acceleration efficiency, magnetization 
level, and coherence length scale all increase in time. 
These trends and the above lower limits on e acc and 
€b indicate that a shock propagating into a cold, 
homogeneous plasma with B = remains a vi- 
able model for astronomical shocks, with no need 
for additional assum ptions about magnetic turbu- 
lence generation (e . g., iGoodman and MacFadvenl 120071 
iMilosavljevic et al.l 120071 : ISironi and Goodman! 12007ft . 
Our results confirm that particle acceleration and mag- 
netization are intimately related, with high energy par- 
ticles playing a major role in generating the magnetic 
fields which in turn scatter and accelerate the particles. 

We thank J. Arons and P. Goldreich for helpful dis- 
cussions. U.K. is a Friends of the Institute for Advanced 
Study member, and is supported by NSF grant PHY- 
0503584. The research of B. K. & E. W. is partially sup- 
ported by AEC, ISF & Minerva grants. A.S. acknowl- 
edges the use of computational resources at TIGRESS 
computing center at Princeton University and the sup- 
port from Alfred P. Sloan Foundation fellowship. 



REFERENCES 



R. Blandford and D. Eichler, Phys. Rep. 154, 1 (1987). 

O. Buneman, in Computer Space Plasma Physics, edited by 

H. Matsumoto and Y. Omura (Terra Scientific, Tokyo, 1993), 

p. 67. 

P. Chang, A. Spitkovsky, and J. Arons, ApJ 674, 378 (2008). 

J. T. Frederiksen, C. B. Hededal, T. Haugbolle, and A. Nordlund, 

ApJ 608, L13 (2004). 
J. Goodman and A. MacFadyen, ArXiv e-prints (2007), 

arXiv:0706.1818. 
A. Gmzinov, ApJ 563, L15 (2001). 

A. Gruzinov and E. Waxman, ApJ 511, 852 (1999). 

C. H. Jaroschek, H. Lesch, and R. A. Treumann, ApJ 618, 822 
(2005). 

T. N. Kato, ApJ 668, 974 (2007). 

B. Katz, U. Keshet, and E. Waxman, ApJ 655, 375 (2007). 

J. G. Kirk, Y. Lyubarsky, and J. Petri (2007), arXiv:astro- 
ph/0703116. 

A. Loeb and E. Waxman, Nature 405, 156 (2000). 

L. Maraschi, in Active Galactic Nuclei: From Central Engine to 
Host Galaxy, edited by S. Collin, F. Combes, and I. Shlosman 
(2003), vol. 290 of Astronomical Society of the Pacific 
Conference Series, p. 275. 



M. V. Medvedev and A. Loeb, ApJ 526, 697 (1999). 

M. Milosavljevic, E. Nakar, and F. Zhang, ArXiv e-prints (2007), 

arXiv:0708.1588. 
T. Piran, Reviews of Modern Physics 76, 1143 (2005). 
L. O. Silva, R. A. Fonseca, J. W. Tonge, J. M. Dawson, W. B. 

Mori, and M. V. Medvedev, ApJ 596, L121 (2003). 
L. Sironi and J. Goodman, ApJ 671, 1858 (2007). 
A. Spitkovsky, in Astrophysical Sources of High Energy Particles 

and Radiation, edited by T. Bulik, B. Rudak, and G. Madejski 

(2005) , vol. 801 of American Institute of Physics Conference 
Series, pp. 345-350. 

A. Spitkovsky, ApJ 673, L39 (2008a). 
A. Spitkovsky, ApJ 682, L5 (2008b). 

E. Waxman, Plasma Physics and Controlled Fusion 48, B137 

(2006) . 



